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Abstract 

The molecular dynamics (MD) simulation study of solvation structure and free energetics in 
l-ethyl-3-methylimidazolium chloride and l-ethyl-3-methylimidazolium hexafluorophosphate using 
a probe solute in the preceding article [Y. Shim, M. Y. Choi and H. J. Kim, J. Chem. Phys. 
xxxxxxxxxx] is extended to investigate dynamic properties of these liquids. Solvent fluctuation 
dynamics near equilibrium are studied via MD and associated time-depedent friction is analyzed via 
the generalized Langevin equation. Nonequilibrium solvent relaxation following an instantaneous 
change in the solute charge distribution and accompanying solvent structure reorganization are 
also investigated. Both equilibrium and nonequilibrium solvation dynamics are characterized by at 
least two vastly different time scales — a subpicosecond inertial regime followed by a slow diffusive 
regime. Solvent regions contributing to the subpicosecond nonequilibrium relaxation are found to 
vary significantly with initial solvation configurations, especially near the solute. If the solvent 
density near the solute is sufficiently high at the outset of the relaxation, subpicosecond dynamics 
are mainly governed by the motions of a few ions close to the solute. By contrast, in the case of a low 
local density, solvent ions located not only close to but also relatively far from the solute participate 
in the subpicosecond relaxation. Despite this difference, linear response holds reasonably well in 
both ionic liquids. 







I. INTRODUCTION 



In the preceding paper,— referred to as I hereafter, we have investigated the equihbrium sol- 
vation structure and free energetics in two different room-temperature ionic hquids (RTILs), 
l-ethyl-3-methyhmidazohum chloride (EMI+Cl^) and l-ethyl-3-methylimidazolium hexaflu- 
orophosphate (EMI"'"PFg") by employing a diatomic and a benzenelike probe solute.— We have 
found that these solvents show a dramatic electrostrictive behavior; i.e., the radial and angu- 
lar distributions of solvent ions around the solute become structured markedly as the solute 
charge separation is increased. This results in, e.g., significant outer-sphere reorganization 
free energies, which are comparable to and sometimes in excess of those in highly dipolar 
solvents.— Also the effective polarity of these RTILs, which measures their solvating power 
of dipolar solutes, is very high. The free energy curves relevant to solvent fluctuations and 
associated local force constants vary with the solute charge distributions. This indicates that 
a linear solvation scheme is generally not valid for RTILs. Nevertheless the harmonic approx- 
imation with state-dependent force constants provides a reasonable framework to quantify, 
e.g., the solvent reorganization free energy. 

In this article, we continue our investigation of solvation in RTILs, focusing on their 
dynamic properties. We study equilibrium solvation dynamics and examine their variations 
with the solute electronic structure. Using the generalized Langevin equation (GLE) descrip- 
tion, we analyze time- dependent solvent friction. As for nonequilibrium, time-dependent 
Stokes shifts subsequent to an instantaneous change in the solute charge distribution and 
accompanying solvent structure relaxation are investigated. Underlying molecular motions 
and their roles in nonequilibrium relaxation are also examined. Comparison is made between 
equilibrium and nonequilibrium solvation dynamics to gain insight into linear response in 
RTILs. 

The outline of this paper is as follows: In Sec. II we briefly review the time-correlation 
function and GLE descriptions of solvation dynamics. The molecular models and simulation 
methods employed in our study are summarized in Sec. III. In Sec. IV we present an analysis 
of the simulation results for equilibrium fluctuation and nonequilibrium relaxation dynamics 
in EMI'^'Cl^ and EMI"'"PFg . Concluding remarks are offered in Sec. V. 
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II. SOLVATION DYNAMICS AND GENERALIZED LANGEVIN EQUATION DE- 
SCRIPTION 



Here we review briefly the time-correlation function description of solvation dynamics 
to establish our notation. We also give a brief account of the GLE equation approach to 
solvation. 

As in I, the solute is assumed to be described by two nonpolarizable electronic states, a 
and h. For a given solvent configuration Q and a solute electronic state i (= a, 6), we denote 
the total energy of the combined solute-solvent system as Ei{Q). The Franck-Condon (FC) 
energy /\Ea^b{Q) associated with the a ^ h transition in the presence of Q is [Eq. (4) of I] 

^E^^,{Q) = E,{Q) - E,{Q) . (1) 

As mentioned in I, AEa^b provides a very convenient variable to describe the collective 
solvent influence on the solute.— From experimental perspective, by measuring the solute 
absorption or fluorescence band, in particular, modulations in AEa^b induced by solvation, 
we can probe solvent structure and dynamics in the presence of the solute in state a. Hence- 
forth, we will refer to the originating state of the probe FC transitions (viz., the state a in 
AEa^b) as the active state, while the other electronic state involved in the FC transitions 
(i.e., b in AEa->b) as the reference state. In the remainder of this paper, we focus mainly 
on dynamics associated with AEa^b under both equilibrium and nonequilibrium solvation 
conditions. 

We first consider the equilibrium solvation situation. Fluctuating dynamics of AEa^b 
near equilibrium in the presence of a solute in state a is usually described by the normalized 
time-correlation function 

{^AE^m^AE^AtTl 

^'^/^^'^ = — mE^^b?) — 

with 

{6AEa^b{<^) 6AEa^b{t)) = J dQo (Qo) 6AEa^b{Qo) 6AEa^b{Qt) , (3) 

where (. . .) denotes an average with the ensemble distribution function /^^ in equilibrium 
with the a-state solute, 6AEa-^b = ^Ea^b ~ (^Ea^b) is the deviation of AEa^b from its 
equilibrium average, and Qo and Qt represent, respectively, solvent configurations at time 
and t. We point out that for the equilibrium ensemble average in Eq. (j21), the solute state 



which determines the equihbrium ensemble distribution at t = and the active state which 
the FC transitions to monitor solvent dynamics originate from are identical. 
The normalized time-dependent Stokes shift function 



AB„i(0) - AB„i(oo) ^ ' 

directly accessible via time-dependent fluorescence,— "^Si is widely used to describe nonequi- 



librium solvation dynamics. In Eq. AEa-*b{t) denotes the average FC energy associated 
with the a ^ b transition at time t after an instantaneous change in the solute electronic 
structure from the charge distribution of the reference state b to that of the active state 
a. The average is taken over the initial distribution of solvent configurations in equilibrium 
with the reference fe-state^ 

A^(t) = J dQo /r(Qo) ^Ea^,{Qt) . (5) 

From Eqs. and (0), equilibrium and nonequilibrium averages satisfy 



AEa^,{0) = -{AE,^a) = -AE,^aioo) . (6) 

We stress that for both nonequilibrium Sa/b(t) and equilibrium Ca/b(t), solvent dynamics 
occur in the presence of the active a-state solute. The major difference between the two 
lies in the character of their ensembles. In the linear response regime, Sa/b{t) and Ca/b{t) 
are characterized by the same dynamic behavior. Its validity in RTILs is investigated in 
Sec. ITVRl below. 

We turn to a stochastic description of solvation dynamics using a generalized Langevin 
equation (GLE). For brevity, the subscripts representing the solute electronic states are 
suppressed in the remainder of this section. The Mori-Zwanzig projection^ on to a set of 
dynamical variables {6AE, 6AE} yields 

6AE{t) = -ujI 5AE{t) - [ dt' ({t - t') 6AE{t') + R{t) , (7) 

Jo 

where the overdot denotes the time derivative, ({t) and R{t) the time-dependent friction and 
random force (scaled by the inertia associated with 6AE motions), respectively, and Ug the 
solvent frequency defined as 

_ mm 



The friction and the random force are related via the so-called second fluctuation-dissipation 
theorem: 

C(*) = . (9) 

({t) describes dissipative AE motions caused by rapidly-fluctuating R{t), while Us charac- 
terizes inertial dynamics of AE'.— It should be noted that the GLE formulation with Eqs. 
©-Q is formally exact. 

It is easy to show that the equilibrium time-correlation function (TCP) in Eq. (P)) 
satisfies^ 

C{t) = -ul C{t) - I dt' Cit - t') C{t') . (10) 







With the Fourier-Laplace transform 

C[z\ = I dt exp[izt] C{t) , (11) 
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we obtain ^ 

C[z]=(-zz + —^^) , (12) 
V -tz + C[z]J 

which shows that C{t) is governed by the solvent frequecy and friction C('^)- For example, 
near t = 0, the friction term in Eq. (fTUj) can be ig nored, so that C{t) ^ -cj^ C{t). This 
yields C{t) ^ exp[— Icj^t^], manifesting that the short-time behavior is governed by inertial 
dynamics. We also consider the TCP 



(^AA(0)^(t)) 



m = , (13) 



where T^f{t) = 5AE{t)+uj'^5AE{t) represents the "non-conservative" force on 5AE{t), viz., 
the sum of frictional and random forces {{J-'j\f 6AE) = 0). The friction is related to C{t) and 
^(t) via 

{u;! -z^)C[z]+tz I 

l + izC\z] \ u^-z' 
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c>i- '"-:^'T7 =ii+^^^Mi *w. (14) 



In Sec. IIVAI below, we use Eq. fll4|l to determine ({t) from the simulation data. 

III. SIMULATION METHODS 

As in I, the simulation cell is comprised of a single rigid solute molecule immersed in 
either EMI+Cl^ or EMI+PPg solvent, consisting of 112 pairs of rigid cations and anions. 



Two different types of solutes, diatomic and benzenelike molecules, are considered. For each 
solute type, we have employed two different charge distributions: For the diatomic solute, a 
neutral pair (NP) with no charges and an ion pair (IP) with unit charge separation (g = ±e) 
are used, where e is the elementary charge. As for the benzenelike solute, we consider an 
axially symmetric charge distribution appropriate for regular ground-state benzene (RB)^ 
and charge-separated dipolar benzene (DB), which is characterized by one electron transfer 
from one carbon to another carbon atom in the para position, compared to RB. Thus the 
dipole moment difference between the RB and DB charge distributions is 13.5 D. For EMI"*", 
we use the AMBER force field^^ for the Lennard- Jones (LJ) parameters and the partial 
charge assignments of Ref. Q. The united atom representation is employed for the CH2 and 
CH3 moieties of the ethyl group as well as for the methyl group (Fig. 1 of I). Experimental 
geometry determined by X-ray diffract ioni^^ is used. For Cl~, we use a = 4.4 A and e/k-B = 
50.4 K with Boltzmann's constant k-Q- PFj^ is described as a united atom with a = 5.6 A 
and e/fcB = 200K. For the details of the solute and solvent model descriptions, the reader is 
referred to I. 

The molecular dynamics (MD) simulations are conducted in the canonical ensemble at 
T = 400 K using the DL_POLY program.— Equilibrium simulations are carried out with 
2 ns equilibration, followed by a 8 ns trajectory from which averages are computed. For each 
reference state, the nonequilibrium quantities are computed by averaging over 400 distinct 
initial configurations generated from its 8 ns equilibrium trajectory and separated by 20 ps in 
time. These initial configurations are used to generate nonequilibrium solvation by changing 
the solute charges from the reference-state to active-state charge distribution instantaneously 
a la Franck- Condon. The subsequent solvent relaxation to new equilibrium is monitored for 
2 ps in each of the initial configurations. 

IV. MD RESULTS AND DISCUSSIONS 
A. Equilibrium solvation dynamics 

The MD results for Ca/b{t) are presented in Fig. Qfor several values of the solvent density 
p. Regardless of the solute and solvent considered, Ca/bif) is characterized by at least two 
different dynamics: ultrafast initial relaxation followed by an extremely slow decay. In 
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the latter was attributed to ion transport via diffusion and the former to small- 
amphtude inertial translations of ions. There it was also pointed out that the long-time 
relaxation is not a single-exponential decay but can be fitted reasonably well with a stretched 
exponential. Despite their relatively high viscosity, short-time dynamics in both EMI+C1~ 
and EMI+PFg fall in a subpicosecond time regime and account for more than 50 % of the 
entire Ca/bif) relaxation. 

We notice in Fig. ^that the early dynamics of Ca/b{t) tend to become faster with the 
growing charge separation (i.e., dipole moment) of the solute active state. For instance, as we 
change the active state charge distribution from NP to IP (and from RB to DB), the solvent 
frequency Ug [Eq. (jS))] which describes the initial relaxation behavior of Ca/b{t), increases by 
by 30-60 %. Dipolar acetonitrile shows a similar trend according to the present and previous 
simulation studies.— The results in Table H] show that two factors contribute to this trend - 
viz., decrease in {{SAEa-,b)^) and increase in {{6AEa-,b)^) ■ In I, the former is ascribed to the 
enhancement in solvent structural rigidity induced by electrostriction as the solute charge 
separation grows. As for the latter behavior, the increase in the electrostatic interaction 
strength, also discussed in I, plays a central role. To see this, we rewrite {{SAEa^bY) as 

{{SAK^,Y) ^ ((5^v„-V,AE,^,)2) = i^(t;„2)((v„AE,^,)2) 

a a 

= E^!^((v»AW) = 5:^((A/rft . (15) 

a a 

where a labels solvent ions, and are the mass and velocity of ion a, and A/™"' 
is the difference in the Coulomb forces on a arising from the a- and 6-state solute charge 
distributions. For both NP IP and IP NP, Af^°'^^ is, up to a sign, just the Coulomb 
force on a arising from the IP charges because there are no electrostatic interactions between 
the solvent and the NP charge distribution. We have, for the sake of simplicity, ignored the 
cation rotations as well as solute translations and rotations in Eq. (fT3j) . The local solvent 
density increase close to the solute with the growing solute charge separation (see Figs. 2-4, 
6 and 7 of I) enhances the solute-solvent electrostatic interactions, viz., the magnitude of 
A/^°"'. This results in the increase of (((5Ai?a_+b)^) according to Eq. ()15|1 . For the RTILs 
and aprotic CH3CN studied here, {{SAEa^b)"^) and {{6AEa^b)^) contribute in a constructive 
way to the solvent frequency trend with the solute active state charge separation. 

Another interesting feature in Fig.^is that Ca/b{t) becomes more oscillatory as the a-state 
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becomes more dipolar. According to our earlier study,— the solvent structure enhancement 
with the increasing solute dipole is responsible for this trend. Briefly, enhanced structure 
yields a tighter potential of mean force for ions, which results in more oscillations with higher 
frequencies, especially for anions (see Fig.|21below). In contrast to the short-time oscillations, 
the long-time relaxation of Ca/b{t) becomes slower as the active-state charge separation 
increases; both the exponent f3 and time constant tq of its stretched exponential fit, i.e., 
exp[— (t/ro)^], decrease (see Table H)). In Ref. 



13l . this was also attributed to the potential 



of mean force, which becomes more rugged with the growing solute charge separation. This 
decelerates ion transport and thus the long-time decay of Ca/hif)- 

In Fig. 121 we decompose AEa^b into the cation and anion components and analyze the 
contributions of their autocorrelations and cross correlation to Ca/b{t)- The normalization of 
these three components are such that the sum of the cation and anion autocorrelations and 
their cross correlation yields Ca/bit). Figs. El^a) andEfb) show that in EMI"'"C1~, short-time 
subpicosecond solvation dynamics are almost completely determined by anions.— Also the 
oscillatory behavior of Ca/b{t) for a = IP arises from Cl~. Because of their small size, Cl~ 
ions can approach rather close to the solute. This makes their Coulombic interactions with 
the solute very strong. For instance, the Coulomb interaction energy between Cl^ and the 
(-I-) site of IP at a separation of 3.5 A [cf. Fig. 2(a) of I] is about — 95kcal/mol. Therefore, 
translational motions of Cl~ located near the solute, even if small in amplitude, are very 
effective in modulating AE, compared with those of EMI"*" with extended structure and 
charge distribution. This also explains why the anion dominance in short-time dynamics is 
reduced considerably in EMI"^PFg (although the anion contributions tend to be still larger 
than the corresponding cation contributions); translations of bulky PFg do not modulate AE 
as much as those of Cl~ because the former are located at a larger distance from the solute 
than the latter. In addition, due to the enhancement in the cation structure around the solute 
in EMI"'"PFg compared to that in EMI"''C1^ [Sec. IIIA of I], EMI+ plays a more significant 
role in the former than in the latter. Another noteworthy aspect is that the relaxation of 
the cross correlation between the anions and cations is slow. This is presumably due to the 
local electroneutrality condition that strongly correlates anions and cations. Finally, we note 
that as in the solvent spectral shifts studied in I (see Table II there), the role of anions in 
equilibrium solvation dynamics becomes somewhat reduced in the presence of benzenelike 
solutes (results are not shown here) compared to diatomic solutes. We nonetheless notice 
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that the general characteristics of Ca/b{t) are very similar between the two solute types. 

The time- dependent friction for solvent dynamics in the GLE description (Sec.|n|) is 
exhibited in Fig. |21 To obtain ({t), we have first computed Ca/b{t) and ^{t) and used their 
Fourier-Laplace transforms in the second and third expressions of Eq. (jl4|) to determine (^[z] 
for small and large z, respectively. These two functions have been connected smoothly, so 
that resulting ([z] when substituted into Eq. ()12|) reproduces C{t) accurately in the entire 
time domain . ^^1^^ The results for EMI''"C1~ in the presence of diatomic solutes show that 
analogous to Ca/b(t), there are two totally different dynamic regimes for ((t), ultrafast relax- 
ation with a time scale shorter than 0.1 ps and an extremely slow decay. Due to numerical 
uncertainties involved in the determination of C{t), it is not easy to determine unambigu- 
ously if the latter is a single-exponential decay. Nevertheless, with a single-exponential 
fit, the characteristic time associated with the long-time tail of ({t) would be ~40ps for 
NP and ~ 100 ps for IP. This reveals that the memory effects on solvation dynamics last 
for an extended period in RTILs. As such, the effective friction on solvent motions varies 
considerably with their time scale and makes the long-time behavior of associated Ca/b{t) 
non-exponential. This could have important consequences for reaction dynamics of various 
charge shift and transfer processes. For example, in the case of electron transfer reactions, 
the effective friction relevant to barrier crossing could vary substantially with the barrier 
frequency in RTILs.— 

To place our results for equilibrium solvation dynamics in perspective, we digress here 
to consider potential limitations imposed by the rigid, united-atom description used in 
the present study. While our comparative analysis in I shows that it captures static sol- 
vation properties well, we expect that solvation dynamics — both ultrafast and long-time 
components — in real ionic liquids would be accelerated compared to our predictions through 
participation of cation (and to some extent, PFq) internal motions that are absent in the 
rigid, united-atom description. To check this, we have conducted a test simulation using 
the flexible, aU-at„.des:.p«„„ of Ref.Qfo.. EMI. a.d*e rigid, all-atom description of 
Ref. Q as in I. In Fig. Ufa), its preliminary results for Ca/b(t), averaged over a 2 ns trajec- 
tory in the presence of the diatomic solute are compared with the results of the united atom 
description. While long-time dynamics of the all-atom description with flexible cations are 
faster than those of the rigid, united-atom model, the presence of clear demarcation between 
the subpicosecond and diffusive regimes is common to both models. Furthermore, the rela- 
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tive contributions (~ 60-70 %) and time scale (~0.5 ps) of the subpicosecond dynamics are 
essentially the same between the two. It is quite remarkable that two models with rather dif- 
ferent parametrization and representation yield this level of quantitative agreement in their 
TCF predictions. One of their most noticeable differences though lies in the initial relaxation 
near t = 0, which is considerably faster for the flexible cation model than for the rigid one. 
For better understanding, we decompose Cip/Np(t) of the former into the cation and anion 
components in Fig. |3fb) just like the rigid case above. Comparison with Fig.|2fd) shows that 
it is indeed the cation component of the flexible model that is responsible for its rapid oscil- 
latory decay in the first ~0.1ps.— The oscillation period is ~30fs, which corresponds to a 
frequency of ~ 1000-1200 cm~^. Its likely sources are bending motions of various CH groups 
of EMI"*". While the internal modes of this frequency may not be excited near equilibrium at 
room temperature, they will nonetheless participate in nonequilibrium relaxation with excess 
energy (see Figs. |H1 and ^1 below). More generally, our classical analysis here, though pre- 
liminary, shows that ion internal motions, e.g., torsion and bending modes of polar groups, 
could be important to energetics and dynamics in RTILs. Thus while the overall aspects 
of Ca/b{t) can be captured quite reasonably by rigid potential models, a proper account of 
ion internal motions will be needed to accurately quantify solvation dynamics and related 
dynamic electronic spectroscopies, e.g., optical Kerr effect spectroscopy.— 

B. Nonequilibrium Solvation Dynamics 

The MD results for Sa/b(t) given by Eqs. (jlj) and (0) are shown in Fig.^l Just like Ca/b(t), 
Sa/b{t) is characterized by at least two distinct relaxation processes in all cases we studied. 
About 50-70% of the entire solvent relaxation occurs in less than 0.3 ps, irrespective of the 
probe solute. This confirms the earlier conjecture based on dynamic and static Stokes shift 
measurement sS-iiO that a significant portion of solvent relaxation occurs on a subpicosecond 
time scale in RTILs. We notice in Fig. El that subpicosecond dynamics of 5'Np/ip(t) (and 
'S'rb/db('^)) are faster than those of 5'ip/Np(t) (and 5'db/rb(''^))- We will return to this point 
later on. 

In Figs. El and 5'^/^(t) and S^^l(t), the anion and cation components of Sa/bit), are 
presented for EMI+Cl^ and EMI+PFg . While both the anions and cations play an important 
role in relaxation dynamics, the contributions from the former are much larger than those 
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from the latter in EMI+Cl". As mentioned above, this is due to the small size of Cl~. For 
EMI"*"PFg with bulky anions, the contributions from the anions and cations are comparable. 
We nonetheless notice that the anion contributions to the subpicosecond solvent dynamics, 
measured as the initial drop in S^J^{t) in the first 0.2 ps or so, tend to be somewhat larger 
than the cation contributions. Even in the case of the benzenelike solutes in Figs. [7|^c) 
and[7|^d) where S^jKO) is larger than S'^J?^(0), the anion contributions to the subpicosecond 
dynamics are still larger than the corresponding cation contributions. Another noticeable 
feature is that as in the Ca/b(t) case, the relative contributions of cations and anions to Sa/b(t) 
vary with the solute type; the cations tend to play a more important role in the presence 
of the benzenelike solute than in the presence of the diatomic solute (see above). In I, it is 
conjectured that variations in solute structure and charge distributions influence respective 
roles of cations and anions and thus details of solvation by modulating the ion distributions 
close to the solute. This generally supports the view that solvation properties of ionic liquids 
measured via dynamic and static electronic spectroscopies can vary with the probe solutes.— 
Nevertheless, overall features of S a/bit) in Fig. El do not seem to be that sensitive to the 
solute types although our investigation is limited to the diatomic and benzenelike solutes. 

Before we proceed further, we briefly pause here to consider MD statistics employed for 
nonequilibrium calculations. We note that at t = 0, the anion and cation components of 
Sa/b{t) satisfy 

S:T(^) = ^,7r°(0) . (16) 
This is due to the fact that despite the nonequilibrium character of Sa/b{t) and its compo- 
nents, their initial values depend only on equilibrium quantities [cf. Eqs. ^ and ©]. We 
notice in Figs. Inland [7| that Eq. is well satisfled by our MD results. This indicates that 
400 conflgurations used to compute averages provide decent statistics for nonequilibrium. 

Returning to our main thread, we decompose cation motions into translation and rotation 
and study their respective contributions to the solvent relaxation . ^^i^^i^^ For this purpose, it 
is convenient to consider AAEa-^b{t) = A.Ea-,b{t) — AEa^bi^), which starts from the null 
value at t = 0, for each type of motion. According to the results in Fig. |H1 the contributions 
from cation translations exceeds those from cation rotations by a factor of two or more, 
revealing that translational motions of ions are of primary importance in S^'^l{t) throughout 
the entire relaxation.— This holds for all cases we have studied. 

In view of huge electrostrictive effects found in I, it is needless to say that the nonequihb- 

11 



rium solvent relaxation is accompanied by significant solvent rearrangements in RTILs. To 
understand the link between the structural relaxation and solvation dynamics, we examine 
the time-evolution of solvent structure subsequent to an instantaneous change in the solute 
charge distribution at t = 0. The results for the evolving radial distributions of the cations 
and anions around the diatomic solute in EMI+Cl^ are exhibited in Fig. El We first consider 
Figs, ini^a) andin^b), where the solvent, initially equilibrated to the IP solute, relaxes to a 
new state in equilibrium with the NP charge distribution. This corresponds to 5'Np/ip(t) 
since the solvent relaxation occurring in the presence of the active NP state is monitored 
via the NP IP electronic transitions. The time evolution of the number of Cl~ in the first 
solvation shell, viz., within 5 A of the positive site of the solute is given in Table ITTl For 
convenience, we will refer to this quantity as the CI" coordination number hereafter. We 
notice that the structural change in the first ~ 0.2 ps is not that significant on a relative 
scale. During this period, the evanescent first solvation shell moves out a little and becomes 
broadened somewhat; its peak height decreases by ~35%. The reduction in the Cl~ coordi- 
nation number is less than unity from 4.1 to 3.4. Only after ~lps, does the destruction of 
the first solvation shell become obvious. Despite this, 5'Np/ip(t) relaxes by more than 60% 
during the first ~ 0.2 ps as shown in Figs. El and In terms of energy, this corresponds to 
a change of ~30kcal/mol in A-E^p^ip [cf- Fig- Eta)]. 

To gain further insight into the ultrafast relaxation of Sa/b{t), we divide anions into the 
first solvation shell and "bulk" region^ at t = and investigate the temporal behavior 
of their respective contributions to AEa-^b- The MD results corresponding to S'Np/ip(t) 
are given in Fig. ITnT a). where the cation results are also presented for comparison. The 
contributions from the first solvation shell anions vary significantly with time, especially at 
small t, whereas those from the bulk change little. This reveals that the contributions to 
subpicosecond dynamics of S^^jpit) arise almost completely from Cl~ in the first solvation 
shell! Though lesser in extent, this is also the case with cations. This paints a picture 
that motions of a few ions close to the solute essentially govern the fast nonequilibrium 
solvation dynamics occurring immediately after the instantaneous IP NP change in the 
solute charge distribution. It is remarkable that these few ions account for energy relaxation 
of ~40 kcal/mol in about 0.2 ps, which corresponds to more than 60 % of the entire 5'Np/ip(t) 
decay. 

We investigate the time evolution of the average distance between the solute and ions. 
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which are initially located in the first solvation shells. The results in Fig. [TTT a) show that 
ion translations are characterized by two different regimes: The displacement at short time 
grows more rapidly than a linear increase in t. We ascribe this to "inertial" translations 
initiated by sudden changes in the Coulomb force on ions exerted by the solute at t = 0. The 
displacement at later time increases nearly linearly with t, suggesting a diffusive transport 
regime for ions, subject to the Coulomb force (thus biased). It should be noticed that the 
average displacement of the first solvation shell ions in the inertial regime is not significant; 
their net displacement during the first 0.2 ps is about 0.5 A. But due to their strong Coulombic 
interactions with the solute, this is sufficient to modulate and relax Ai?NP/ip by more than 
60 %. We thus conclude that in the case of the initial IP NP change in the solute electronic 
state and ensuing solvent relaxation, rapid inertial translations of small amplitude associated 
with the ions situated very close to the solute are responsible for the ultrafast subpicosecond 
dynamics of 5'Np/ip(t).^^ The subsequent slow decay of 5'Np/ip(t) is attributed to (biased) 
diffusive migrations of these ions into the "bulk" region^ as well as diffusive rearrangements 
of bulk ions. 

We now turn to the opposite case, where the solvent, initially equilibrated to the NP 
solute, readjusts to the the newly-created IP charge distribution. Analogous to the 5'Np/ip(t) 
case analyzed above, the initial relaxation of the Stokes shift function 5'ip/Np(t) is extremely 
fast. The t-dependent solvent distributions in Figs. IHfc) and|nfd), on the other hand, show 
that the formation of solvent structures around the solute is quite sluggish compared to the 
short-time relaxation of S'ip/Np(t). To be specific, the nascent first solvation shells develop a 
descernible structure only after ~0.4ps. During this period, the Cl^ coordination number 
changes just a little from 1.8 to 2.15 (see Table |n}. In the first ~ 0.2 ps, during which 
more than 50% of the Sip/^p(t) relaxation is completed [Fig. Efa)], the increase in the 
Cl~ coordination number is mere 0.2! For comparison, recall that the Cl~ coordination 
number change during the same 0.2 ps period is ~0.7 in the case of S'Np/ip(t). Thus while 
rapid inertial translations of ions (see Fig. fTlT b)) are expected to be responsible for the 
subpicosecond relaxation of 5'ip/Np(t) as in S'Np/ip(t), the ions close to the solute alone may 
not be sufficient to account for the ultrafast energy relaxation of 30-40 kcal/mol. 

To pursue the last point above a little further, we group the anions together in three 
different regions according to their separation from the solute at t = and analyze their 
contributions to AEa^b{t) in Fig. llOf b). For convenience, we refer to the three regions as 
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"near" (r < 5 A),— "intermediate" (SA < r < TA), and "far" {r > 7 A) regions, where r is 
the initial distance between the solute positive site and Cl~. For completeness, the cation 
results are also shown there.— We notice that anions from all three regions make substantial 
contributions to subpicosecond solvation dynamics. For example, the contribution from Cl~ 
in the far region exceeds 10% of the total anion contribution to ultrafast >S'ip/np(^)- This is 
in contrast with the 5'Np/ip(t) case above, where short-time dynamics are almost completely 
governed by ion motions in the near region [Fig. ITUT a)]. This is due to the difference in 
the anion density in the near region between the two cases (cf. Fig. In the 5'Np/ip(t) 
case, the number of Cl~ initially in the near region is high, so that the large magnitude 
of their Coulombic interactions with the solute dominates Aii^NP^ip- By contrast, in the 
case of 5'ip/Np(t), due to the absence of Coulombic interactions between the NP solute and 
ions, there are on the average fewer than two Cl~ ions present near the solute at the outset 
of the solvent relaxation. These are not sufficient to dissipate the extra Coulomb energy 
generated by the change in the solute charge distribution at t = 0. Therefore the Cl~ in the 
intermediate and far regions play a considerably more signficant role in the subpicosecond 
relaxation of S'ip/Np(t) than in the case of 5'Np/ip(t). This also explains why the former are 
somewhat slower than the latter as we noticed in Fig. above. Since the electrostatic force 
on ions arising from the solute decreases rapidly with their separation, it will take longer for 
ions at a large distance from the solute to dissipate energy through translations than those 
close to the solute. Slow relaxation of S^p^j^p(t) occuring on a time scale longer than 1 ps 
mainly involves migration and redistribution of Cl^ originating from the intermediate and 
far regions. 

To summarize, while the basic features of the Stokes shift functions 5'Np/ip(t) and 5'ip/Np(t) 
are very similar, the range of Coulombic interactions that determine the energy gap relaxation 
differs considerably between the two. To be specific, in the case of 5'Np/ip(t), rapid receding of 
the first solvation shell ions from the solute essentially determines the subpicosecond solvation 
dynamics; thus the solute-solvent Coulombic interactions at close ranges are of primary 
importance. By contrast, in the 5'ip/np(^) case. Fig. fTUr b) shows that ions at all distances 
contribute to the solvation dynamics through their initial approach to the solute ( "clustering 
in" ) via small amplitude inertial translations. Therefore the Coulombic interactions between 
the solute and solvent relevant for S'ip/Np(t) dynamics are of considerably longer range than 
those for 5'NP/ip(t). (For an analysis of related issues in dipolar liquids, see Ref.y.) While 
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our conclusion is based exclusively on the analysis of EMI+Cl", it also applies to EMI+PFg . 
The major difference between the two is the relative importance of cations; i.e., they play a 
more significant role in the latter than in the former. 

Before we conclude, we consider the validity of linear response in RTlLs. Comparison 
between Ca/b(t) and Sa/bif) in Fig- IT^ shows that despite the large electrostrictive effects 
observed in I (Figs. 2-6 there) and the difference in the range of Coulombic interactions 
relevant to Sa/b{t) found above, linear response holds surprisingly well in EMl+Cl^. In 
particular, Cip/Np(t) and Sjsip/ip(t) yield an excellent agreement in short-time dynamics and 
similarly for other pairs, Cdb/rb(^) and S'rb/db(^), etc. This is presumably due to similarity 
in the solvent configurations Qt between Ca/b(t) and Sb/a(t) for small t [cf. Eqs. OJandE], so 
that their early dynamics occur in a similar solvent environment. (Recall that the two are 
characterized by the same ensemble distribution function, i.e., their solvent configurations 
are the same at t = 0.) In a similar vein, we would expect that at long time, slow dynamics of 
Sb/aif) would resemble those of Cb/ait) (rather than those of Ca/bif)) because Qt associated 
with Sb/aif) eventually evolves to the configurations equilibrated with the solute 6-state (see 
Fig.iniabove). While this appears to be the case in Fig. El we need to perform nonequilibrium 
simulations considerably longer than 2 ps to test this notion unambiguously. Finally, we note 
that the electronic polarizability absent in the current study — especially that of the probe 
solutes — will introduce a significant nonlinear behavior in solvent response.— 



V. CONCLUDING REMARKS 



In the present article, we have studied dynamic properties of solvation in EMI^'^Cl^ and 
EMI"'"PFg via MD simulations. We have investigated dynamics Stokes shifts subsequent to 
an instantaneous change in the solute charge distribution and accompanying solvent struc- 
ture reorganization that evolves in time. We have studied fluctuating solvent dynamics near 
equilibrium and analyzed associated friction dynamics via the generalized Langevin equa- 
tion. It was found that both equilibrium and nonequilibrium dynamics show two totally 
different relaxation processes: rapid subpicosecond dynamics arising mainly from inertial 
ion translations and slow relaxation ascribed to ion transport. While the solvent dynamics 
are dominated by Cl~ motions in EMI+Cl^, both anions and cations play an important role 
in EMI+PF^. We have also found that linear response holds reasonably well in both ionic 
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liquids. 

We have examined molecular details of nonequilibrium solvent relaxation by analyzing 
cation and anion motions and their contributions from different solvation regions. We found 
that cation translations are considerably more important than cation rotations during the 
entire solvent relaxation.— The range of solvation regions which participate in the subpi- 
cosecond relaxation varies significantly with the initial solvent density near the solute. In 
the case of a high local density near the solute at t = 0, ensuing subpicosecond solvent re- 
laxation is mainly governed by the motions of a few ions close to the solute. By contrast, in 
the opposite case of a low initial density, solvent ions not only close to but also relatively far 
from the solute contribute to the subpicosecond relaxation dynamics. Despite this difference, 
their dynamic Stokes shift functions show very similar characteristics. 

It would be worthwhile in the future to extend the present study based on a rigid, united- 
atom description in more realistic directions. Among others, the incorporation of molecular 
flexibility into simulations in the all-atom representaion would be desirable. Our preliminary 
results considered here indicate that though overall dynamic features are well captured by 
the rigid model, it clearly misses several important details, in particular, at short time, 
due to the absence of ion internal motions, such as torsion and bending. Another aspect 
not included in our model description is electronic polarizability of ions. While its effect 
in RTILs may not be as significant as in polar solvents, its proper account is needed to 
accurately quantify e.g., solvent spectral shifts and reorganization free energy. 

Finally, solute dynamic properties, such as rotational friction and vibration energy re- 
laxation, in RTILs and their variations with solute electronic structure will be reported 
elsewhere.— 
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In the case of optical probe molecules, AE^™/^ + Ai?^" represents the Stokes shift associated 
with the electronic transitions between a and b. 

Within the numerical accuracy of 1 kcal/mol which is the bin size employed for AEa^h to 
calculate Pa/b and Fa/b (Eqs. (??) and (??)), our MD results for the A^;^"^, values and {AEa^b) 
are the same. 

In view of the discussions above, it would be the opposite in the case of water, e.g., Ajp/^p is 
smaller than Ajyfp/ip. 

The first solvation shells for anions and cations (and their corresponding bulk regions) are defined 
with respect to the positive and negative site of the dipolar solute, respectively. 
These ions mainly correspond to the anions and cations of the first solvation shells, which are 
defined with respect to the positive and negative sites of the solute. 
^'^ By decomposing the cation motions into rotational and translational components, we have in- 
vestigated their respective contributions to AE relaxation. We have found that the contribution 
from rotational motions is considerably larger than that from translations. 
However, there is a caveat: the internal motions of anions and cations, e.g., bending and torsional 
motions of charge and polar groups of EMI"*" , absent in the current rigid model description can 
play an important role, especially in very early dynamics [cf. Fig. 

the internal motions of cations are totally absent in the rigid model description we employed. In 
Sec. IIVAI above, we observed that the initial relaxation of Ca/bit) [Fig. E] of the flexible model 
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and r > 8 A, respectively, where r is the distance between the cation center and solute (+) site. 
The bulk is defined to be the rest of the anions (cations) not included in the first solvation shell. 
Though neglected here, anion internal motions, i.e., PF stretch and bend, would introduce 
similar rapid oscillations (of ~ 40-60 fs) to the anion component and thus to Cip/np(^)- We 
nonetheless believe that their contribution would be smaller than that of cation internal motions 
because for a given frequency, the average amplitude of internal motions associated with more 
massive F is considerably smaller than that associated with lighter H. 
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TABLE I: MD results'^' 



solvent 


density 


a/6 




{i6AEa-.br) 


(PS~^) 


/?'') 




EMi+cr 


1.1 


NP/IP 


83.3 


7520 


9.5 


0.25 


0.31 






IP/NP 


60.5 


14720 


15.6 


0.11 


0.15 


EMi+cr 


1.1 


RB/DB 


77.5 


7140 


9.6 


0.5 


0.76 






DB/RB 


54.2 


13530 


15.8 


0.1 


0.49 


EMI+Cl- 


1.2 


NP/IP 


72.3 


10590 


12.1 


0.12 


0.14 






IP/NP 


57.9 


15010 


16.1 


0.05 


0.005 


EMI+PFg 


1.31 


NP/IP 


79.8 


2870 


6.0 


0.25 


0.69 






IP/NP 


64.5 


3820 


7.7 


0.15 


0.13 


EMI+PFg 


1.375 


NP/IP 


77.8 


3290 


6.5 


0.2 


0.54 






IP/NP 


52.0 


4310 


9.1 


0.1 


0.077 


EMI+PFg 


1.375 


RB/DB 


81.2 


4810 


7.7 


0.6 


2.22 






DB/RB 


38.7 


5120 


11.5 


0.05 


0.0935 


CHgCN^) 


0.73 


NP/IP 


46.4 


4640 


10.0 










IP/NP 


41.3 


6870 


12.9 







"^Units for density, AEa^b and time (t and Tq) are g/cm^, kcal/mol and ps, respectively. 
^^A stretched exponential fit to Ca/bit), i-e., exp[— (t/ro)'^], for t < lOOps after ultrafast 
initial relaxation. 

^^MD results at T = 300 K. The LJ parameters and partial charges employed are from Ref. 0. 
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TABLE II: Time evolution of the CI eoordination number around the diatomic solute after the 
change in the solute charge distribution 



active state 


tips) 
site 


0.0 


0.1 


0.2 


0.4 


0.6 


0.8 


1.0 


2.0 


oo 


NP 


(+) 


4.07 


3.41 


3.40 


3.35 


3.28 


3.19 


3.04 


2.55 


1.80 




(-) 


0.34 


0.43 


0.71 


1.25 


1.58 


1.54 


1.71 


1.79 


1.80 


IP 


(+) 


1.80 


1.77 


2.00 


2.15 


2.31 


2.39 


2.46 


2.87 


4.07 




(-) 


1.80 


1.45 


1.22 


1.02 


0.88 


0.86 


0.86 


0.65 


0.34 
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FIG. 1: Equilibrium time correlation function Ca/b{t) ii^ (a) EMI^Cl at density p = l.lg/cm^ 
and in (b) EMI^PFg at /? = 1.375 g/cm^ for the solute active/reference states a/b: NP/IP ( — ); 
IP/NP ( ); RB/DB (•••); DB/RB ( ). 

FIG. 2: Components of Ca/bit) for (a) NP/IP and (b) IP/NP active/reference states of the diatomic 
solute in EMI+Cl" and for (c) NP/IP and (d) IP/NP in EMI+PF^. The dotted and dashed hues 
represent the contributions from the cation and anion autocorrelations, respectively, while the 
dashed-dot line denotes the negative of their cross correlation divided by 2. The RTIL densities 
are the same as those in Fig. ^ 

FIG. 3: Time-dependent solvent friction (^(t) (in units of ps~^) associated with Ca/bit) in EMI+Cl" 
at p = 1.1 g/cm^: a/b = NP/IP ( — ) and IP/NP (•••). The corresponding initial values C(0) are 
marked with x and □, respectively. Details of the initial behavior are shown in the inset. 

FIG. 4: (a) Comparison of Ca/b{t) between the united-atom (UA) and all-atom (AA) descriptions for 
EMI+PFg . The solvent densities are p = 1.375 g/cm'^ (UA) and 1.34 g/cm'^ (AA). (b) Components 
of Cip/Np(t) in the all-atom representation in (a). 

FIG. 5: Stokes shift function Sa/b{t) in (a) EMI+Cl" and (b) EMI+PF^. The solute ac- 
tive/reference states employed are NP/IP (— ), IP/NP ( ), RB/DB (• • •) and DB/RB ( ). 

The densities are the same as those in Fig. ^ 

FIG. 6: Components of Sa/b{t) in EMI+Cl^ for the solute active/reference states a/b: (a) NP/IP; 
(b) IP/NP; (c) RB/DB; (d) DB/RB. 

FIG. 7: Components of S^/bit) in EMI+PFg' for the solute active/reference states a/b: (a) NP/IP; 
(b) IP/NP; (c) RB/DB; (d) DB/RB. 

FIG. 8: Contributions from rotational and translational motions of ions to the relaxation of 
AEa^b{t) (in units of kcal/mol) after the initial (a) IP NP and (b) NP IP changes in the 
solute charge distribution in EMI+PFg . The corresponding dynamic Stokes shift functions for (a) 
and (b) are 5'NP/ip(t) and 5'ip/np(0i respectively. 
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FIG. 9: Time evolution of ion radial distributions subsequent to an instantaneous change in the 
solute charge distribution in EMI+Cl^. The anion and cation distributions given, respectively, in 
(a) and (b) describe the structural relaxation accompanying S^p/ip{t), while the corresponding 
distributions in the opposite Sip/-!<jp{t) case are shown in (c) and (d). For clarity, plots of the radial 
distribution at later time are shifted upwards from the initial distribution according to its time 
evolution in each pane and only the results for the methyl group of EMI"*" (i.e.. Ml in Fig. 1 of I) 
are presented in (b) and (d). 

FIG. 10: Contributions to AEa^b{t) (in units of kcal/mol) from different solvent regions for 
(a) NP/IP and (b) IP/NP active/reference states of the solute. The solid and dotted lines repre- 
sent the anion and cation contributions, respectively. In (a), EMI+ ions located within 6 A from 
the negative site of the solute are defined as the first solvation shell, while the 5 A separation is 
employed to define the corresponding first solvation shell for Cl~ around the solute positive site. 

FIG. 11: Time evolution of the average distance r (in units of A) of the ions from the solute after 
the instantaneous (a) IP NP and (b) NP — > IP changes in the solute charge distribution. Only 
the ions located initially in the first solvation shells are considered. 

FIG. 12: Comparison of equilibrium Ca/bit) and nonequilibrium Sa/b{i) in EMI^Cl" in the presence 
of (a) diatomic and (b) benzenelike solutes. 
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